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The effect of asymmetry in tlie transbilayer lipid distribution on the dynamics of phase separation 
in fluid vesicles is investigated numerically for the first time. This asymmetry is shown to set a 
spontaneous curvature for the domains that alter the morphology and dynamics considerably. For 
moderate tension, the domains are capped and the spontaneous curvature leads to anomalously slow 
dynamics, as compared to the case of symmetric bilayers. In contrast, in the limiting cases of high 
and low tensions, the dynamics proceeds towards full phase separation. 
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Asymmetric distribution of lipids in the two leaflets 
of the plasma membrane is ubiquitous to many eu- 
karyotic cells. Most of phosphatidylserine and phos- 
phatidylethanolamine are located in the cytoplasmic 
leaflet, while sphingomyeline and phosphatidylcholine are 
predominantly in the outer leaflet jl|. This asymme- 
try, maintained by the cell through many active and 
passive processes, plays an important role in the lat- 
eral and trans-membrane compositional and morpho- 
logical organizations in the nanometer scale. A very 
good example for such organization is the nanoscale do- 
mains, referred to as rafts. These domains are believed 
to be liquid-ordered regions, mainly composed of sphin- 
gomeylin, which is a saturated lipid, and cholesterol [3. 
In spite of the wealth of experimental studies on lipid 
rafts, the mechanisms leading to their formation and 
their stability remain under intense debate. These issues 
are complicated by the presence of many components and 
processes in biomembranes. 

With the aim to achieve the understanding of the phys- 
ical properties of biomembranes, many experimental |3| 
and theoretical |l,|a] investigations have been carried out 
on relatively simple model lipid membranes. To the best 
of our knowledge, in all these studies, membranes have 
the same lipid composition in both leaflets of the bilayer. 
The natural next step in complexity towards the under- 
standing of biomembranes is to consider membranes with 
different lipid composition in the two leaflets. An impor- 
tant question that arises is then: What role does this 
transbilayer asymmetry play on the domain structure of 
lipid bilayers? In particular, will this asymmetry result 
in a flnite size of these domains? The purpose of the 
present paper is to address these questions using large 
scale dissipative particle dynamics (DPD) simulations. 

Recent experiments clearly demonstrated that in mul- 
ticomponent lipid vesicles there is a strong registration of 
the domains in the two leaflets. That is the type of lipids 
in the two leaflets are locally the same. This implies that, 
if there is a transbilayer asymmetry in the average lipid 
composition, these domains in register in the two leaflets 
will have to be of different areas. In order to minimize 
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FIG. 1: (a) The configuration that minimizes the interaction 
energy between the lipids of a bilayer when the domains in the 
two layers are of the same area and (b) when the areas of the 
domains are different. The spontaneous curvature resulting 
from this area difference can be estimated by Eq. (1). 



the interaction energy between unlike lipids the domains 
have to curve as shown schematically in Fig. 1. The 
area difference between the domains in the two leaflets 
should then give rise to a spontaneous curvature. This 
spontaneous curvature is only a function of the average 
compositions in the two leaflets, (jjout and (f>in defined as 
the number of A- lipids in the outer and inner leaflets, 
respectively, divided by the total number of lipids in the 
same leaflet 
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where, e is the bilayer thickness. 

Mean fleld calculations of multicomponent monolayer 
membranes, with curvature coupled to the local concen- 
tration through a spontaneous curvature, predict that 
when one of the components has a preferred curvature, 
the equilibrium state is that of caps or stripes. These 
calculations, considering an ordered patterns of domains, 
have been carried out both in two- 0, H IE S U^ and 
three-component '8| monolayer models for lipid mem- 
branes. An ordered array of buds has also been proposed 
in the strong segregation limit |3, HOl ■ These calculations 
are performed on an infinite planar membrane. An ex- 
tension of these calculations to the case of finite closed 
vesicles are difficult and have not yet been reported. It 
is important to note that the Hamiltonian used in these 



monolayer models are based on a local bilinear coupling 
of the form, J. (j)H, where H is the local mean curva- 
ture and (f) is the local composition field. One way to 
obtain these models from two-component bilayer models 
amounts to having a local mismatch between types of 
lipids in the outer and inner leaflets 0)13' ^^^ therefore 
absence of domain register. Instead, a strong register re- 
ported in experiments suggests that it is more reasonable 
to expect a spontaneous curvature generated by the dif- 
ference in area between the domains that are in register. 

Motivated by the argument above, we carried out a 
systematic DPD simulation, of self-assembled lipid bi- 
layer vesicles with different lipid compositions in the two 
leaflets of the bilayer. Within the DPD approach [Ij| a 
number of atoms or molecules are coarse-grained in order 
to form a fluid element, thereafter called a dpd particle. 
Speciflcally, we have in our case three types of beads, cor- 
responding to a water-like bead, labeled w, a hydrophilic 
bead labeled h, representing a lipid head group, and a hy- 
drophobic bead, labeled t. representing a group of CH2's 
in the tail group jlSl [lq |. The model parameters are 
selected such that the membrane is impermeable to the 
solvent thus allowing, in the case of a vesicle, to inves- 
tigate the effect of conservation of inner volume. As in 
our previous study a lipid particle is modeled by a fully 
flexible linear amphiphilic chain, constructed from one 
hydrophilic /i-particle, connected to three consecutive hy- 
drophobic i-particles 4, jll] . There are two types of lipid 
particles, corresponding to A and _B-lipids. The heads 
and tail dpd particles of an a-lipid, with a = A or B, are 
labeled ha and ta, respectively. The position and veloc- 
ity of each dpd-particle are denoted by r j and v^ , respec- 
tively. Within the DPD approach, all particles are soft 
beads that interact with each other through three pair- 
wise forces corresponding to a conservative force, pf^', a 
random force, f}"', and a dissipative force, f.<°'. The con- 
servative force between dpd particles, i and j is given by 
^If '= aijuj{rij)nij where Vij = r^ - Tj and n,y = ^ij/rij. 
Since all dpd particles are soft, we choose w(r) = 1 — r/rc 
for r < re, with Vc is the cutoff of interaction and is 
used to set a length scale in the simulations. w(r) = 
for r > Vc- The integrity of a lipid chain is ensured 
through a simple harmonic interaction between consecu- 
tive monomers, F^^'= —C{l — \rij\/b)nij, where the spring 
constant C = lOOe and the length b = 0.45rc. The 
dissipative force, describing friction between neighbor- 
ing particles, is given by F(f)= -7w^(|rij|)(ny • Vij)ny, 
where v^ — v.^ — Vj . Finally, the random force is given 
by, F(f>= -crw(|ry|)Cjj(At)~i/2ny where Cij is a ran- 
dom noise with zero mean and unit variance, and At 
is the time step of the simulation. The fluctuation- 
dissipation theorem requires that 7^ = afj2k-QT. We 
used a fixed value for the parameter Cy — a for all pairs 
of dpd-particles. The head-head and tail-tail interac- 
tions between unlike lipids are given by ahA.hs — 50e 
and atA,tB = 50e respectively with e setting a scale for 
energy. All other parameters used in the present study 
are the same as in references \4. Illj. In our simulations. 




FIG. 2: Snapshot sequence of System II (bottom row) as com- 
pared to that of a symmetric vesicle having the same area to 
volume ratio (top row). Snapshots from left to right corre- 
spond to i = 100, 2000 and 4000r, respectively. 



a vesicle is composed of 16 000 lipid particles, embedded 
in a fluid consisting of 1 472 000 solvent particles with 
density p — 3r^^, and the time step. At = O.OSr [I3 
with T — {rarll eY'"^ and ra is the mass of a single dpd 
particle. 

The interactions above ensure that the membrane is 
impermeable and that flip-flop events are extremely rare. 
The tail-tail interaction between two like lipids is less re- 
pulsive than that between two unlike lipids. Different 
tail-tail interactions are needed in order to mimic the fact 
that the hydrophobic region in domains which are rich in 
the saturated lipid and cholesterol is conformationally 
different from domains rich in unsaturated lipids. The 
phase separation is initiated through randomly labeling 
a fraction of all lipid particles as types A or B, such that 
the outer and inner leaffets have compositions (pout and 
(pin, respectively. The coarsening dynamics is then mon- 
itored through the computation of the cluster sizes of 
the minority lipid and the net interfacial length between 
the segregated domains. By simultaneously monitoring 
both interfacial length and the average domain area, we 
are able to investigate the physical mechanism via which 
domain growth proceeds. 

Here, we present results for the case of (<?!>o«t,0m) — 
(0.4, 0.2) with three different area-to-volume ratios, cor- 
responding to the case of the number of solvent par- 
ticles in the vesicle's core, Ns — 138000 (System I), 
Ns = 112400 (System II), and N, = 72300 (System III). 
In the remaining of this letter, the symmetric case, refers 
to a vesicle with {4>out,4>in) = (0.3,0.3) and an area- 
to-volume ratio equal to that in System II. In Fig. 2, 
a sequence of snapshots in the case of system II, with 
{Pout, Pin) — (0.4,0.2) are compared to that of a vesicle 
with symmetric transbilayer composition corresponding 
to {pout. Pin) = (0.3,0.3). This flgure clearly shows that 
the onset of domain capping is shifted to earlier times by 
the transbilayer asymmetry in the lipid composition, as 
opposed to the case of a symmetric vesicle. This capping 
result from the register between domains in the outer 
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FIG. 3: Net interface length as a function of time. The data 
with the slope —4/9 corresponds to System II and that with 
slope -0.13 corresponds to the symmetric vesicle. Inset shows 
the net number of domains as a function of time for System 
II (top curve at late times) and for the symmetric vesicle 
(bottom curve). 



and inner leaflets and mismatch between their areas (see 
Eq. CD)- 

At late times the coarsening dynamics is slower in the 
asymmetric case than in the symmetric case as is evi- 
dent from Fig. 2. This is substantiated by the time de- 
pendence of the net interface lengths of the two systems 
shown in Fig. 3. This slowing down is not observed in 
the symmetric case even during late times, which the 
expected fast growth law, L ^ t^'*/^, resulting from co- 
alescence of caps 3. Note that while the curvature of 
domains in the symmetric case is set by the competi- 
tion between line tension and bending modulus 5], in 
the asymmetric case it is induced by the spontaneous 
curvature resulting from the asymmetry in the compo- 
sitions of the two leaflets, and is set at early times (see 
Fig. 2). The late times slowing down shown in Fig. 3, 
must be the result of this spontaneous curvature. It is 
interesting to note that during the late stages, the interfa- 
cial length has a very small growth exponent ('^ —0.13), 
indicative of a logarithmic growth. This is normally at- 
tributed to microphase separation, ubiquitous to many 
other soft materials ,17|. This non-algebraic slow dynam- 
ics is the result of the competition between interfacial 
tension, which is the driving force of the phase separa- 
tion, and an effective long-range repulsive interaction. In 
the present case this repulsive interaction should result 
from the combined effect of spontaneous curvature and 
lateral tension. Microphase separation in multicompo- 
nent lipid bilayers with spontaneous curvature has been 
predicted by mean field theories both in the weak and 




• •* 



FIG. 4: Snapshot sequences for the three area to volume ratios 
considered in the present study. Rows from top to bottom 
represent the Systems-I, II, and III, respectively. Snapshots 
from left to right, in each row, correspond to times i = 100, 
1000, 2000, and 5000r, respectively. 




FIG. 5: Net interface length versus time. Curves from bottom 
to top corresponds to Systems I (bottom). III (middle), and II 
(top), respectively. The data line for System I (bottom) hass 
shifted downwards for clarity. Straight lines are guides to the 
eye. In the inset, the average interface length per domain is 
shown versus time in Systems II (top curve), and III (bottom 
curve) . 



strong segreg ation limits at low values of lateral tension 
[y,0,|a,l3ill3- We should note that these theories assume 
that domains, in the microphase separation regime, are 
organized into a regular lattice. However thermal fluctu- 
ations could easily destroy the long-range order of these 
domains. Indeed, in our simulation, we see a melt of caps 
instead of a regular lattice. 

Surface tension is expected to play an important role 



on the phase separation of multicomponent lipid mem- 
branes. In the case of closed and impermeable vesicles 
with conserved number of lipids, the role of surface ten- 
sion is effectively played by the area-to- volume ratio. Ex- 
perimentally, surface tension is controlled by an osmotic 
pressure between the inner and outer environments of 
the vesicle. In Fig. 4, snapshots of vesicles with asym- 
metric composition are shown for three different values 
of the area-to- volume ratio, corresponding to Systems I, 
II and III. At low tension (System III), the large amount 
of excess area allows domains to acquire a bud shape, 
such that interfaces of the corresponding domains in each 
leaflet are in register. Due to similar reasons, at interme- 
diate tensions (System II) domains acquire a cap shape 
since the excess area available in this case is less than 
that in System III. At high tensions (System I), the lack 
of excess area prevents domains from curving, leading 
them to remain nearly flat throughout the phase separa- 
tion process. 

The net interface lengths of the domain structure ver- 
sus time are shown in Fig. 5 for the three different ten- 
sions. This figure shows that during early times of the 
phase separation process, i.e. for t < 250 the dynamics 
is independent of tension. As demonstrated by the early 
time configurations in Fig. 4, this is due to the fact that 
during this regime, domains are flat, and the composition 
dynamics is decoupled from the membrane curvature dy- 
namics. At late stages, the dynamics in the three systems 
become noticeably different. We notice in particular that 
the dynamics is very slow in the case of intermediate ten- 
sion (System II), as discussed above. In contrast, domain 
coarsening is faster and algebraic for both high (System 
I) and low (System III) tension, implying an approach to- 
wards full phase separation. A growth law, L{t) ~ t~^/^ 
in System I, is the result of coalescence of flat domains 
as is shown earlier in Refs. [Jillll- Our results for cases 
of intermediate and high tensions are in agreement with 
recent theoretical predictions JJlJ. However, for the low 
tension case, we see growth dynamics towards full phase 
separation. The equilibrium state here will be a com- 
pletely phase separated vesicle 18]. 



The diffusion coefficient of a single bud should scale as 
D ~ 1/a^/^, where a is the area of a single bud. If domain 
coarsening at late stages of system III is mediated by coa- 
lescence of spherical buds, with a fixed neck radius Z, dif- 
fusing on the surface of the vesicle (the length scale here 
being set by the ratio of difference in area occupied by the 
B phase to the bilayer thickness) then the net interface 
length, L{t) = N{t)l -- t~^/^, where N is the total num- 
ber of buds on the vesicle ^^ . This growth law is much 
faster than what is shown in Fig. 5. l{t) = L{t)/N for 
Systems II and III plotted in the inset of Fig. 5. This fig- 
ure shows that the interfacial length per bud, and hence 
the neck radius, in system III reaches saturation at about 
t « 250t. a close examination of snapshots at late times 
clearly indicate that domain growth in system III is medi- 
ated by coalescence. Therefore, in order to reproduce the 
observed time dependence, L{t) '-^ t'^'^/'^ , the buds must 
move with an effective diffusion coefficient that scales as 
1/a^. Since this is unrealistic there must be another pro- 
cess opposing the necks of two buds from merging. The 
details of this process is currently under investigation. 

In conclusion, we examine the effect of transbilayer 
asymmetry in the lipid composition on phase separation. 
We found that at intermediate tension, the asymmetry 
leads to anomalously slow coarsening of caps, in agree- 
ment with recent mean field calculations. At low or van- 
ishingiy small lateral tension, where the domain structure 
is that of buds, and at high tension, where the domains 
are fiat, we found algebraic domain growth leading to full 
phase separation. It would be very useful to explicitly 
correlate the nature of domain growth, with asymmetric 
transbilayer lipid distribution, with the surface tension 
of the membrane. We plan to perform the calculations of 
the surface tension in the various regimes from the lateral 
stress profiles using the approach similar to that in Refs. 
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